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We describe a simple model of evolution which incorporates the branching and extinction of 
species lines, and also includes abiotic influences. A first principles approach is taken in which 
the probability for speciation and extinction are defined purely in terms of the fitness landscapes 
of each species. Numerical simulations show that the total diversity fluctuates around a natural 
system size JVnat which only weakly depends upon the number of connections per species. This is in 
agreement with known data for real multispecies communities. The numerical results are confirmed 
by approximate mean field analysis. 
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The Bak-Sneppen model was introduced to illustrate 
the possible role of self-organised criticality in evolving 
ecosystems (Bak 1993, 1997). It is a toy model that de- 
scribes every species by a single scalar quantity, relating 
to the expected time before that species evolves to a new 
form. Interactions in the ecosystem are incorporated by 
assuming that a species that evolves can alter the time 
taken for other species to evolve, such as those involved in 
predator-prey or host-parasite relationships. The model 
is said to be self- organised critical because it approaches 
a critical state without any apparent need for fine pa- 
rameter tuning. Consequently, it predicts that extinction 
events of magnitude AE should occur with a frequency 
proportional to (AE)~ a , which is consistent with known 
paleobiological data (Sole 1996). Further evidence has 
come from analysis of the temporal distribution of ex- 
tinctions, which exhibits "1/f noise", also predicted by 
the model (Sole 1997a). 

Other simple models of macroevolution have now been 
devised which also claim agreement with the paleobio- 
logical data (Peliti 1997). Some of these are believed 
to be self-organised critical (Sole 1997b), although oth- 
ers exhibit power law behaviour via different mechanisms 
(Roberts 1996, Newman 1997). All of these models have 
in common the assumption of a constant system size. 
This has been justified by assuming that each species oc- 
cupies a single ecological niche, and that if a species is 
made extinct its niche is immediately filled by a simi- 
lar, newly emerged species. The concept of an ecolog- 
ical niche refers to a set of conditions and interactions 
within the ecosystem that only a single species can sat- 
isfy. However, since the definition of a niche depends 
upon the other species, the total number of niches should 



be defined from within the system itself rather than be- 
ing fixed to some arbitrary constant value for the benefit 
of computer simulation. 

Models have been devised which are similar to the Bak- 
Sneppen model but allow the total number of species to 
vary in time. Kramer et al. introduced a model in which 
the species are placed onto a branching phylogenetic tree 
structure, where each species only interacts with its clos- 
est relatives (Kramer 1996). Depending upon the choice 
of a parameter, the tree either expands indefinitely or 
stops growing after a finite time. In the model of Wilke 
et al. the system refills after an extinction event at a 
rate according to a parameter iV max , which is "the max- 
imal number of species that can be sustained with the 
available resources" (Wilke 1997). However, since the re- 
sources are themselves biotic, we believe that any such 
N max should be defined from within the system. 

In this paper, we derive and study a version of the Bak- 
Sneppen model in which the probabilities for speciation 
and extinction are defined purely in terms of the individ- 
ual species. Nonetheless, the total diversity fluctuates 
around a natural system size without the need for global 
control. In particular, there is no recourse to "ecologi- 
cal niches" . The rules of the model are based on careful 
considerations of the motion of species on their fitness 
landscapes. This is described in Sec. [n|, and the results 
of numerical simulations of the model are presented in 
Sec. Jnj. Comparisons to real ecosystems are made in 
Sec. Approximate analysis of the model is presented 
in Sec. [V] which supports and enhances the numerical 
work. Finally, we discuss our results in Sec. VI. 
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II. QUANTITATIVE DESCRIPTIONS OF 
MACROEVOLUTION 

To quantitatively describe an evolving ecosystem re- 
quires some general principle that applies equally to all 
species. A candidate for such a principle is to consider 
the relationship between an organism's genotype and its 
fitness, which is some measure of the expected number of 
genes passed back into the species' gene pool (Dawkins 
1983). Each point in the multidimensional space of all 
possible genotypes can be assigned its own fitness value, 
forming a fitness landscape, as schematised in Fig. [j]. 
The process of evolution can then be described as a walk 
over this landscape in the direction of increasing fitness. 
Rather than try to calculate the fitness for each genotype 
from first principles, clearly a hopeless task, Kauffman 
has assumed that the relationship is so complex that it 
can be well approximated by random variables (Kauff- 
man 1993). This leads to the concept of rugged fitness 
landscapes, where "rugged" refers to the large variations 
in fitness that can result from small changes in genotype. 
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FIG. 1. Schematic of a typical fitness landscape for a 
species labelled i, where for clarity the space of all possible 
genotypes has been compressed onto a single axis. The fitness 
scale is arbitrary. The barriers have heights of fn and fa. 

Models based on the Bak-Sneppen approach assume 
that each species moves on its own rugged fitness land- 
scape, eventually becoming trapped in the region of a lo- 
cal maximum. If the landscape is fixed, then the species 
can only evolve by moving to a different maximum. Sup- 
pose that a species labelled i is at a local maximum 
which is separated from nearby maxima by fitness bar- 
riers of heights , where j = 1,2,3,... and the fo are 
ordered such that fo < fo+i- Over time, fluctuations 
in the species' fitness may bring it into the vicinity of 
one of its barriers, allowing it to crossover to a differ- 
ent maximum. This constitutes an evolutionary event 
in which the species changes from one typical genotype 
to another. For uncorrelated fluctuations, the expected 
time Tij to crossover a barrier of height fo will be given 
by an Arrhenius equation of the form 

Tij ~ exp(/y// ) , (1) 



where the constant fo fixes the timescale and is analogous 
to temperature. 

Since the landscape is rugged, a species that escapes 
from one maximum will soon become trapped by another, 
and will find itself surrounded by an entirely new set of 
barriers fo. In the limit fo — > 0, ([!]) implies that it will 
always be the species i with the smallest fo in the sys- 
tem that evolves first, and that the other species will have 
moved no appreciable distance towards their own barriers 
by the time this occurs. Thus the ecosystem will consist 
of species that infrequently hop between maxima at a rate 
that depends upon the fo , but are otherwise essentially 
static. The evolution of species i will alter the landscapes 
of all those species k linked to it in the ecosystem, for in- 
stance via predator- prey or host-parasite relationships. 
Although each species k will in general have to move on 
their new landscapes before finding a new maximum, it 
is a further approximation of the theory that this can be 
ignored and only the barriers fkj are affected. 

The original Bak-Sneppen model is defined purely in 
terms of the smallest barriers fn , which are arranged on 
a lattice in such a way that interacting species occupy ad- 
jacent lattice sites. The evolutionary process described 
in the previous paragraphs then reduces to the dynami- 
cal interaction between adjacent fn. The system is static 
until the site with the smallest fn evolves, when fn and 
all of the ff-i in adjacent sites k are assigned new val- 
ues. The system is again static until another site evolves, 
and so on. It has been shown that the essential system 
behaviour is insensitive to details such as the choice of 
probability distribution for the fn (Bak 1993, Paczuski 
1996). This robustness relates to the universality of the 
critical state, and is essential if such a simple model is to 
faithfully describe real ecosystems. 

The Bak-Sneppen model can be enhanced by a more 
detailed consideration of the fitness landscapes and their 
interaction. Three features absent from the original 
model will be considered here, namely speciation, extinc- 
tion and external noise. Each feature is described in gen- 
eral terms below before the new model is fully specified. 

Speciation: Speciation occurs when two sub-populations 
reach a state of reproductive isolation and should be con- 
sidered as separate species (Maynard-Smith 1993, Ridley 
1993). For instance, two groups that are reunited af- 
ter prolonged geographical isolation may have evolved so 
much in different directions that they are unable to pro- 
duce viable offspring. Up until now, a species has been 
described as simply occupying a region of genotype space. 
More detailed analysis shows that a population forms a 
"cloud" of points of roughly equal fitness around the local 
maximum (Kauffman 1993). Normally the whole popu- 
lation crosses over the same barrier fn, but if fo ~ fn 
then it is possible that part of the population will in- 
stead cross over the barrier fo. If this happens, the two 
subpopulations will move to different maxima and a spe- 
ciation event will have occurred. A simple criterion for 
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speciation is to say that species i will branch into two 
subspecies if 

fa - fn < Ss (2) 

when it evolves to a new form, where 5s is a constant 
parameter. Further barriers could also be considered to 
incorporate the simultaneous splitting into three or more 
subspecies, but such events will be very rare and are ig- 
nored here. 

Extinction: The system size would increase without limit 
if only speciation were allowed, so some mechanism is re- 
quired by which a species can be made extinct and per- 
manently removed from the system. The original Bak- 
Sneppen model does not distinguish between this form 
of extinction and pseudo-extinction, which is where a 
rapidly evolving species disappears from the fossil record 
if its intermediate forms are not recorded. What is re- 
quired is some criterion for true extinction defined purely 
in terms of the individual species' fitness landscapes, 
analogous to (§) . It is not clear how this may be achieved. 
Instead, a heuristic approach is adopted here, which is to 
say that a species k is made extinct if it is linked to the 
species with the minimum barrier, and has fki greater 
than some threshold value. This proves to be the sim- 
plest choice for which the system size does not diverge. 

External noise: A fitness landscape is ultimately a func- 
tion of the species itself, the species with which it in- 
teracts, and any factors external to the ecosystem, so 
fluctuations in the inorganic environment can also cause 
the fitness barriers to change. Examples include local 
disturbances such as volcanic eruptions or the formation 
of a new river, to global events including meteor impacts 
and changes in the sea level. The interactions between 
species have already been incorporated into the model, 
but no allowance has yet been made for these abiotic 
factors. Continuing with the philosophy that changes in 
fitness can be approximated by random variables, exter- 
nal influences are assumed to alter the fitness landscapes 
by an amount O(Sg) per unit time, where 5g is a new 
parameter. More precisely, every species in the system 
will have their barriers altered by an amount 

fij -> fij + Sgij , (3) 

where the Sgij(t) are uniformly distributed on [—4^, 4f] 
and uncorrelated in time. External effects will occur on 
a separate timescale to the evolutionary processes in (Q) , 
but for simplicity both timescales are fixed at the same 
constant rate in this model. 

It remains to be decided how interacting species are 
linked together. The original Bak-Sneppen model placed 
the species on a regular crystalline lattice in which inter- 
acting species occupy adjacent sites, but this is not flexi- 
ble enough to incorporate the addition of new species to 



the system and is of no use here. Real food webs have 
a much more involved structure, and if the full range of 
interactions is allowed rather than just links in the food 
chain, then it appears that a great many species interact 
at least weakly (Hall 1993, Caldarelli 1998). Trying to 
model this would only serve to draw attention away from 
the main motivation for the new model, which is to allow 
a variable system size. Instead, we adopt the mean field 
approach in which each species i interacts with K — 1 
other species k chosen at random from the system. The 
species k are reselected at every time step. 

The extended model can now be fully specified. 
The ecosystem consists of N(t) species labelled by 
i = 1 . . . N(t). Each species occupies the region around 
a local maximum on a rugged fitness landscape, and is 
separated from nearby maxima by barriers of various 
heights f^. The larger barriers can be ignored since 
they will rarely contribute to the dynamics, but at least 
two must be retained if speciation is to be incorporated. 
Hence each species is defined by its two smallest barri- 
ers fi\ and fi2, which are uniformly distributed over the 
range [0,1] and then ordered so that fa > fn. The fol- 
lowing steps (i)-(vi) are iterated for every time step. 

(i) Evolution: The smallest fn in the system is found and 
marked for evolution. It will move to a new maximum in 
step (vi). 

(ii) Speciation: If the single species marked for evolution 
has fi2 — fn < Ss , then a new species is introduced to 
the system with random barriers. N — ► N + 1. 

(iii) Interaction: K — 1 other species are chosen at ran- 
dom from the remaining N — 1 in the system. They will 
be assigned new barriers in step (vi). 

(iv) Extinction: If any the of the K — 1 interacting species 
has fn > 1, it is removed from the system. N — > N — 1. 

(v) External noise: Every barrier /y in the system is 
transformed according to (||) , and reordered if necessary. 

(vi) New barriers: The species marked for evolution in 
step (i) and the K — 1 interacting species from step 
(iii) are assigned new random barriers, ordered such that 
fi2 > fn- 

With such specific definitions of general processes, it 
is obviously important to check that the model is robust 
to any arbitrary choices. To test this, the simulations 
have been repeated with various changes to the rules, 
and in no case was any qualitative deviation observed. 
For instance, different values for the extinction threshold 
in step (ii) give the same behaviour, even if the threshold 
value varies in time around a fixed mean. Both 5s and Sg 
were chosen from uniform, Gaussian and exponential dis- 
tributions, again with no apparent change in behaviour. 
Since the model appears to be robust, further discus- 
sion will be restricted to the simple set of rules given in 
(i)-(vi). The threshold value for extinction was fixed at 
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1 to minimise the number of new parameters. 

Before continuing, it should be pointed out that the 
algorithm presented in steps (i)-(vi) above is not exactly 
the same as that described in our previous exposition 
of this work (Head 1997). This earlier model assumed 
that all K species "mutate" (evolve) at every time step, 
whereas it is of course just the species with the mini- 
mum barrier that evolves. The corrected model stud- 
ied here behaves in much the same way as its previous 
incarnation, except that the number of species is now 
only weakly dependent on the connectivity K. This is in 
agreement with data for real multispecies communities, 
as discussed further in the next section. 



III. RESULTS OF NUMERICAL SIMULATIONS 

The quantity of interest is the total system size N(t). 
This varies in a manner that depends upon the choice 
of values for the parameters K, Ss and Sg, as described 
below. 

• Ss = Sg = 0: Steps (ii), (iv) and (v) never feature in 
the time evolution of the system and the larger bar- 
riers fi2 are redundant. The fn interact in the same 
way as the original Bak-Sneppen model, the only differ- 
ence being that each fn is the smaller of two uniform 
distributions on [0,1] and so is distributed according to 
P(fa) = 2(1 - f a ), f a e [0, 1]. Since the model is ro- 
bust to the choice of probability distribution, this differ- 
ence is not important. N(t) remains fixed at its initial 
value. 

• K = 1: There arc no interactions, and the species that 
evolves will always have fn < 1 so extinction is impossi- 
ble. N(t) — > oo if Ss > or remains fixed if Ss = 0. 

• K > 1, Ss > and Sg — 0: There is speciation but no 
extinction, N(t) — > oo. 

• K > 1, Ss = and Sg > 0: There is extinction but no 
speciation, N(t) —> 0. 

• K > 1, Ss > and Sg > 0: N(t) fluctuates around 
some constant value iV nat which is independent of the 
initial system size. An example is given in Fig. |^. Note 
that if Sg ^> Ss, N nat is so small that statistical fluctu- 
ations will eventually send N(t) — > and the simulation 
is over. 

That -/V na t should exist at all is by no means obvious, 
since N(t) does not explicitly appear in the rules for spe- 
ciation and extinction. It exists because of the external 
noise of order Sg, which is just as likely to push two barri- 
ers apart as to bring them together and so does not affect 
the rate of speciation. However, the noise acts asymmet- 
rically on barriers near the threshold for extinction, tend- 
ing to push species over this threshold into the small tail 
corresponding to those species that will be made extinct 
when next selected. Since every species is subjected to 



external noise at every time step, the rate of extinction 
increases with N whilst the rate of speciation remains 
roughly constant. A steady state will be found when 
these two rates balance. This qualitative reasoning is 
confirmed by the analysis in Sec. M. 
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FIG. 2. Plot of N(t) against t starting from a system with 
200 species. K = 4, Ss = 0.008 and Sg = 0.02. 



IV. COMPARISON TO REAL MULTISPECIES 
COMMUNITIES 

The parameters Ss and Sg are abstract quantities de- 
fined purely in terms of the model, so it is not possible 
to estimate their values for real ecosystems. Nonetheless 
it is intuitively reasonable to assume that speciation and 
extinction events are rare, and therefore both Ss and Sg 
should be small. The number of links per species K — 1 
has been measured for real communities, and according 
to some studies is independent of the system size (Hall 
1993, Kauffman 1993). This is in agreement with numer- 
ical simulations of the model, which shows only a weak 
dependence on K from the range K = 2 to K = 16, as 
given in Table | (at end of document). There is a small 
peak around K ss 4, which also corresponds to the most 
common value of K observed in nature. However, the 
data for real ecosystems is based on food webs whereas 
the Bak-Sneppen approach considers all direct interac- 
tions between species, so it is not clear how far this com- 
parison can be taken. 

Turning to consider global ecosystems, the fossil record 
for all marine organisms highlights the possibility of a sta- 
tistical steady state throughout much of the Palaeozoic 
era. The total number of (families of) species fluctuates 
around a roughly constant value up until the mass extinc- 
tion at the end of the Permian period, after which the 
system size increased beyond its earlier levels and is still 
increasing today (Benton 1995). It could be conjectured 
that the new species that emerged after the end-Permian 
extinction were on average either more likely to speci- 
ate, or less susceptible to external noise, or both, which 
should result in an increased system size according to the 



4 



model. The data for continental organisms is less clear 
and if anything shows a continuing increase in diversity 
at varying rates. 

The distribution of the magnitude of the change in 
N per unit time is Poisson to first order, implying that 
the speciation and extinction events are uncorrelated for 
N(t) N nat . However, the distribution is not precisely 
Poisson, which is presumably due to the tendency for 
N to drift towards N nat when! N — N nat | is large. An 
example is presented in Fig. 0. That the distribution 
was not power law is disappointing, but perhaps unsur- 
prising given that the interactions between species are 
randomised at every time step, making it difficult for the 
system to self-organise. It is possible that a spatially 
extended model might allow for correlations to build up 
towards a critical state and a power law to be recovered, 
but this must remain as speculation at present. 




Jump size I N(t+At)-N(t) I 

FIG. 3. Distribution of the absolute change in system size 
per At = 10 3 time steps for K = 2, 8s = 0.008 and Sg = 0.04. 
4 x 10 4 points were sampled over 4 separate runs. 



V. ANALYTICAL DERIVATION OF AWt 

It is possible to derive the dependence of N Dat on the 
parameters K, Ss and Sg by extending the mean field 
solution of the original Bak-Sneppen model (Flyvberg 
1993). In theory this approach could give the exact solu- 
tion, since the interacting species are selected at random 
in the new model and so it is mean field by definition. 
However, the increased dynamical complexity means that 
only the first order parameter dependence has be calcu- 
lated. 

Standard model with two barriers per species 

The original solution was based on a single barrier per 
species. Before tackling speciation and extinction, it 
must first be shown how the mean field approach can be 
modified to handle pairs of barriers. Define p(x, y) dxdy 
to be the probability that a randomly selected species 
has one barrier in the range [x,x+dx) and the other in 
[y,y+dy). Note that this refers to the barriers /y before 



ordering, so x can be less that or greater than y. The 
probability for a species to have both barriers greater 
than to = min(x,?/) is represented by Q(x,y), which is 
related to p(x, y) by 



Q(x,y) 



p(x',y') dx'dy' 



(4) 



Since p(x,y) — for values of x or y outside the range 
[0,1], Q(x, y) = 1 for to < and Q(x, y) = for to > 1. 
The species with the smallest barrier can be any of the 
N in the system, as long as all of the remaining N — 1 
species have larger barriers. Hence p m i n (x,y), the prob- 
ability distribution for the species with the smallest bar- 
rier, is given by 



Pmin[x,y) = Np(x,y)Q N 1 (x,y) . 



(5) 



At each time step, p(x, y) will change by an amount 
Ap(x,y) which is given by the master equation 



&P(%,V) = -JrPmin(x,y) 



K-l ( . . 1 . A K 
- w — l [p{x,y)--p m U^y))+^ 



(6) 



The first term on the right hand side of accounts for 
the evolution of the species with the lowest barrier, the 
second for the new barriers assigned to the K — 1 species 
with which it interacts, and the third term handles the 
K new pairs of barriers. In the statistical steady state, 
Ap = and, using (pi) and M) , 
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N 



K - 1 
N - 1 



P 



N -K 
N - 1 



■pQ 



N-l 



= 0. 



(7) 
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The solution to ([?]) depends upon the behaviour of Q 1 
asiV^oo (Flyvberg 1993). If Q < 1 - 0(1/N), then 
the term proportional to Q N ~ X vanishes and 



p(x,y) 



K 



K - 1 



O 



(8) 



Conversely, if either x or y is so small that p(x, y) = 
0(1/ N), then the second term in (0) will be 0(1/ N 2 ) 
and Q(x,y) = 1 + 0(1/N), so 



PQ 



and hence from 



Pn 



N-l 



K 

N 
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N 2 



K + O 



N 



(9) 



(10) 



Each solution applies in different regions of the (x,y) 
plane, which, for large N, will be separated by sharply 
defined boundaries. These boundaries can be found by 
remembering that both p and p min are probability distri- 
butions and normalise to one. To first order in 1/N, the 
full solutions are 
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p(x,y) 



x an d y > 1 
otherwise, 





A" 



x and y > 1 
otherwise. 



K-l 

k ■ 



K-l 

K 1 



(11) 

(12) 



Hence the species with the smallest barrier will always be 
found in the region where p m i n {x, y) ~ A, and its K — 1 
interacting species will always come from the region cor- 
responding to p(x, y) « Kj (A — 1). 

Analysis for 5s and Sg non-zero but small 

When either 5s > or Sg > 0, the system size N be- 
comes a function of time and the algebra quickly becomes 
prohibitive. The simpler and more intuitive approach 
adopted here is to initially ignore speciation and extinc- 
tion altogether and only incorporate the external noise 
of order 5g. This leads to new solutions for p(x,y) and 
Pmin{x, y), from which the rates of speciation and extinc- 
tion can be calculated even though they arc no longer 
dynamically involved. The natural system size -/V nat is 
then the value of N for which the two rates balance. The 
analysis presented below assumes that 5g is small; large 
values of 5g and 5s are considered at the end of this sec- 
tion. 

The effect of the external noise will be to perturb the 
solutions for p and p min given in JTT| ) and ( |l2| ) , as schema- 
tised in Fig. The master equation (^) must be modified 
in two ways. First, the external noise can cause barriers 
to move outside of the range [0,1], so the range of pos- 
sible x and y must be extended. However, the barriers 
are still assigned values in the range [0,1] and the term 
for new barriers must be altered accordingly. Secondly, 
a term for the noise itself must be included. The new 
steady state equation is 



1 



^e(x)0(l-x)9(y)0(l~y)-^- Y 



P 



N -K 
N — 1 



pQ 



N-l 



Sg 2 i 



(13) 



where the integral is over the region p^ in Fig. |. Strictly 
speaking, the distribution in this equation should be 
p — jjPmin, but this distinction can be ignored for large 
N. When both barriers are large, Q N_1 ~ and ([IT 
can be simplified by the transformation 



(16) 
(17) 



x - 




= a{\ — x) , 


y - 


-y' 


= "(i - y) , 

K - 1 


p- 


^p 


= K P < 




a 2 


48(K-l) 




5g 2 N 



to give 



2 v' 2 PV, y') = P'(x',y') - 6{x ')%') 



(18) 
(19) 

(20) 



For either x' or y' negative, corresponding to x > 1 or 
y > 1, the second term on the right-hand side of (|2(i| ) 
vanishes and the equation can be solved by separation 
of variables. Coupled with the boundary conditions 
p'(x' , y') — > for x' — > — oo or y' — > — oo, the solution is 



p'(x',y') 



ce W+y') 



(21) 



where c is an arbitrary constant. Whatever the value of 
c, it must be independent of K and 5g since these pa- 
rameters do not appear in (p0[). Transforming back into 
the original variables gives the explicit parameter depen- 
dence, 



P(x,y) 



K 



K - 1 



- 2 ( x +v 2 ) for x > 1 and y > 1 . 



Substituting this into (|15|) gives 

A' 



k E oc 5g z N 



A- 1 



(22) 



(23) 



where 6{x) | ^ otherwise. 

The theta functions in the first term ensure that new 
barriers lie in the range [0,1]. The last term on the right 
hand side of ( |l3| ) accounts for the external noise, where 
V 2 is the Laplacian operator. A full derivation of this 
term is given in the appendix. 

Rate of extinction: Each of the A — 1 random neighbours 
will be made extinct if it has x > 1 and y > 1. Thus the 
rate of extinction fee is given by 

/oo poo 
J p(x,y)dxdy, (15) 



Rate of speciation: It has not been possible to find a solu- 
tion to (|20|) for x < 1 and y < 1. The variable separable 
solution does not behave correctly, and other methods 
tried have been fruitless. Instead, the Sg = solution 
will be used as a first approximation. The rate of specia- 
tion ks will be proportional to the density of species with 
\x — y < 5s. Since only the species with the minimum 
barrier can speciate, 

fcs = yy 0{Ss-\x~y\) Pl[Lin (x,y)dxdy . (24) 

Substituting the explicit expression for p m i n (_x,y) from 
(H) into (g|) gives 

k s w 5s (25) 



G 



for small Ss. With Sg > 0, p m i a (x, y) broadens and so 
the real rate of speciation will decrease for larger Sg. 

The value of N na _ t can now be found up to parame- 
ter dependence. The rates of speciation and extinction 
balance when Ue = ks, and therefore 



Nnat OC 



5s_ 
6? 



(26) 



This implies that N nat Sg 2 / Ss should be roughly constant. 
This quantity has been calculate d fr om the numerical 
simulations and is shown in Table III. The agreement is 
good for variations in Sg, but less so for K and Ss. This 
most probably reflects the first order approximation used 
in deriving ks in (p5|). 



1-V 



K-1 




FIG. 4. Graphical representation of the effects of external 
noise on the distribution of barriers p(x,y). Denser shading 
corresponds to higher values of p(x,y). The region labelled 
"Pe" refers to species with both x > 1 and y > 1, which are 
liable to extinction. 



Either Ss or Sg large 

For the sake of completeness, the equivalent expression 
to (26) will now be derived for large Ss or Sg, although 
such values bear no relevance to actual systems. If Ss 
is large but Sg small, the system size rapidly increases 
and with it the expected time a species will move un- 
der the influence of external noise before being assigned 
new barriers. Similarly, if both Ss and Sg are large, then 
the system behaviour is also dominated by the external 
noise. This is called the noise dominated regime. If Ss 
is small but Sg large, N nat becomes so low that fluctua- 
tions will eventually drive every species in the system to 
extinction. 

In the noise dominated regime, p(x, y) will no longer 
be just a small perturbation around the original solution 
but will extend to large positive and negative values in 
both the x and y directions. Since the external noise is 
isotropic, p(x, y) will be symmetric about the x and y 
axes and at most I in 4 species have both barriers in the 



Pe region. Hence the rate of extinction will approach its 
upper bound value of 



K - 1 



(27) 



When a barrier is assigned a new value in the range [0,1], 
it undergoes an unbiased random walk until it is again 
assigned a new value and brought back to near the origin. 
The average number of steps in this walk will be 0(N/K) 
and, since the average step size is 0(Sg), an analogy with 
a one-dimensional random walker implies that the total 
distance travelled will be 0{5g^/N/K) (Papoulis 1991). 
This gives the width of the barrier distribution in both 
the x and y directions. The number of species in the in- 
finite strip x — y \ < Ss is inversely proportional to the 
width of p(x, y), so the rate of speciation is now given by 



(28) 



As N increases, the rate of extinction will remain roughly 
constant but now the rate of speciation will decrease un- 
til a balance is found at N — N nat . From ( p7|) and (|2~F 
the corresponding value of N is 




1 f Ss 



(29) 



A convenient way to display the crossover in behaviour 
from small Sg to the noise dominated regime is to con- 
sider iV na t as a function of S = Ss = Sg. According to (|2^ ) 
and (p9[), this should change from N ~ S" 1 for small S to 
N ~ <5" for large S. Numerical results in support of this 
prediction are presented in Fig. ||. 




In(5) 

FIG. 5. Plot of iVnat versus 5 — Ss — Sg for K = 2, demon- 
strating the crossover to the noise dominated regime. Each 
point corresponds to a single run of 10 6 time steps. The 
dashed lines have slopes of —1 and 0. 
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VI. DISCUSSION 

In summary, we have postulated one possible way 
in which models of macroevolution based on the Bak- 
Sneppen approach can be extended to incorporate speci- 
ation, extinction and abiotic influences. The speciation 
and extinction mechanisms are defined purely in terms 
of each individual species' fitness landscape, irrespective 
of the total number of species in the system. Nonethe- 
less, the total diversity fluctuates around a constant value 
N na t , which was termed the natural system size to stress 
that it was not arbitrarily chosen. 

Although the proposed mechanism for speciation, ie. 
a population simultaneously crossing two different fitness 
barriers, seems appealing, the extinction mechanism is 
far more heuristic and somewhat unsatisfactory. A bet- 
ter model might focus on trying to find a more plausi- 
ble means of extinction, defined in terms of the fitness 
landscapes. For instance, the species chosen for evolu- 
tion might be made extinct if its fitness barrier is below 
some threshold value. It may also be possible to place 
the model on a web structure, and allow the connections 
themselves to be subject to alteration whenever a species 
evolves to a new form. 

The value of iV na t was found to be only weakly depen- 
dent upon the average number of connections per species 
in the system, in agreement with known data. This leads 
us to hope that simple models such as ours may be able 
to reproduce the essential behaviour of real ecosystems. 
More realistic models consider the full fitness landscapes 
rather than just the barriers, but the increased complex- 
ity limits the system sizes that can be simulated (Kauff- 
man 1993). A practical step forward might be to reduce 
known biological principles to simple rules that may be 
applied to global ecosystems. 
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In this appendix, the term for external noise that ap- 
pears in ( |l3| ) is derived. Assuming that Sg is small, noise 
effects alone will result in p(x,y) taking the mean value 
of the surrounding square with sides Sg, that is 



+ 



>ise p(x,y) = 

2 r x+Sg/2 



Sg 2 



x-Sg/2 Jy-Sg/2 



p(x,y) 

y+Sg/2 



p{u, v) dudv . 



(30) 



Since 5g is small, p{x, y) can be expanded according to 
Taylor's theorem as 
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p(x + Sx,y + Sy) = p(x, y) + Sx^ + Sy^- 



Sx 2 d 2 p 
~2\~dx 2 ~ 



SxSy 



d 2 p 
dxdy 



Sy 2 d 2 p 



(31) 



On substituting this into (J30|), the terms in Sx, Sy and 
SxSy integrate to zero, leaving just the leading-order term 

A noisc p(x,y) = ^- ^ 2 p(x,y) + 0(Sg 3 ) , (32) 

where V 2 is the two-dimensional Laplacian operator, 



V 



dx 2 



dy 2 



(33) 



The new term A no j se p(x, y) is added to (j^), the expres- 
sion for Ap(x, y) for when Sg = 0, to give the total change 
in p(x, y) at every timestep for Sg > 0. Setting this to 
zero then gives the new steady state equation (|13|). 
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0.016 
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1257(25) 
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TABLE I. Observed values of the natural system size iV na t 
for different K, 5s and Sg. The standard deviation of the 
fluctuations are given in brackets, which also serve as rough 
error bars. The data has been averaged over at least three 
separate rune of 10 6 — 10 7 timesteps each. Note that the line 
for K = 4, Ss = 0.008 and Sg = 0.02 appears three times to 
allow for easier comparison. 







